Coexistence and critical behaviour in a lattice model of competing species 
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In the present paper we study a lattice model of two species competing for the same resources. 
Monte Carlo simulations for d =1, 2, and 3 show that when resources are easily available both 
species coexist. However, when the supply of resources is on an intermediate level, the species with 
slower metabolism becomes extinct. On the other hand, when resources are scarce it is the species 
with faster metabolism that becomes extinct. The range of coexistence of the two species increases 
with dimension. We suggest that our model might describe some aspects of the competition between 
normal and tumor cells. With such an interpretation, examples of tumor remission, recurrence and of 
different morphologies are presented. In the d = 1 and d — 2 models, we analyse the nature of phase 
transitions: they are either discontinuous or belong to the directed-percolation universality class, 
and in some cases they have an active subcritical phase. In the d — 2 case, one of the transitions 
seems to be characterized by critical exponents different than directed-percolation ones, but this 
transition could be also weakly discontinuous. In the d = 3 version, Monte Carlo simulations are in 
a good agreement with the solution of the mean-field approximation. This approximation predicts 
that oscillatory behaviour occurs in the present model, but only for d > 2. For d > 2, a steady state 
depends on the initial configuration in some cases. 
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I. INTRODUCTION 

Competition is a fundamental force shaping almost ev- 
ery biosystem It operates at the level of species, and 
leads to development of new adaptations and creation or 
extinction of species. The competition occurs also among 
individuals of the same species, e.g., when two neighbour- 
ing plants have to share the same resources like water or 
nutrients in the soil. Also at the cellular level the com- 
petition is present, but various mechanisms usually keep 
it under control. These mechanisms sometimes, however, 
fail and a group of tumor cells with abnormal reproduc- 
tion and differentiation pattern emerges. The growth of 
tumor is a very complex process since tumor cells com- 
pete with normal cells for food, space or waste removal, 
modify vascular system or other tissues, and often lead to 
the death of an organism [2|, |3| ■ Therefore, realistic mod- 
els of tumors should take into account numerous factors. 
Due to this complexity, the multiscale modeling is very 
often used [1] , but even relatively simple models that use 
generalized Lotka-Volterra equations [5| are difficult to 
analyse and understand. 

Since tumor development can be regarded as a spatio- 
temporal pattern-formation process, cellular automata 
seem to provide a natural platform to model such phe- 
nomenon This approach was used, for example, 
to study general aspects of Gompertzian tumor-growth 
problem [13| . or more realistic analysis of the three- 
dimensional brain-tumor model on Voronoi tessallation 
[lij . Some hybrid models that combine cellular automata 
with partial differential equations were also used to de- 
scribe interactions between a tumor and the immune sys- 
tem of the host organism 15| , or morphologies of an avas- 
cular tumor [l6j]. One should also mention that related 
mathematical and computational approaches have been 



used to model of some properties of tissues or genomes 
(SHU- 

The cellular automata implemented in such models are 
also very complex, often with non-local, heterogeneous 
or state-dependent transition rules. It is thus very diffi- 
cult to understand their behaviour on general grounds. 
And some basic insight into the behaviour of such sys- 
tems would be indeed desirable since it could help us 
to understand the role of for example fluctuations and 
competition in such systems. Fluctuations and compe- 
tition play major role in determining the symmetry and 
the nature of the ordering in many statistical mechanics 
systems. They are also known to have spectacular conse- 
quences in chemical or biological systems (2l| . To explain 
phenomena like spontaneous cancer remission, that puz- 
zles medicine for decades [22j , the coupling of fluctuations 
and competition certainly has to be understood. 

Cellular automata can also be successfully applied in 
ecology, where competition and fluctuations are of great 
importance. One of the classical problems in this field is 
related to the competitive exclusion principle formulated 
fifty years ago by Gause (23|. It states that when dif- 
ferent species compete for the same resources, only the 
fittest one survives. However, this principle seems to be 
in contradiction with the observed abundance of species 
coexisting in ecosystems, of which plankton is a partic- 
ularly evident example 24]. Some resolutions of this 
so-called plankton paradox were proposed referring to 
spatio-temporal inhomogeneities [25[, cyclic dominance 
26], or size-selective grazing [13] ■ Taking into account 
the basic nature and potential importance of Gause's 
principle and a variety of interactions in ecosystems, it 
would be desirable to examine this problem also in some 
other models. 

In the present paper we examine a simple lattice model 
of two species competing for renewable resources. The 
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model can be regarded as a generalization of a certain 
lattice prey-predator model 28]. Such models can be 
classified as some kind of cellular automata, but due to 
their asynchronous dynamics they are also closely related 
to some other statistical-mechanics models, namely par- 
ticle systems |29|. The species in our model obey the 
intraspecific (but not interspecific) exclusion rule. Us- 
ing Monte Carlo simulations and a mean-field approxi- 
mation, we examine phase diagrams and phase transi- 
tions, which mark the limits of existence or coexistence of 
species. The Monte Carlo simulations are made for one-, 
two-, and three-dimensional versions of the model, while 
the mean-field approximation enabled us to make some 
predictions concerning the coexistence of species or the 
oscillatory behaviour even for higher-dimensional mod- 
els. The two species in our model can also be interpreted 
as normal and tumor cells and actually, we will use such 
a terminology throughout our paper. In this context we 
show examples of tumor remission and recurrence, and 
of different morphologies of invading tumors. One has 
to be aware, however, that such an interpretation must 
be taken with considerable care since the present model 
provides only a very crude description of competition be- 
tween normal and tumor cells. 



II. MODEL 

Lattice models are frequently used to describe various 
problems in ecology [30] . Such individual-based approach 
very often supplements more traditional techniques based 
on differential equations [3lJ • A classical problem in this 
field, having origin in the works of Lotka and Volterra, are 
prey-predator systems [H,|33, 36]. An intensive research, 
especially in the physicists community, is inspired by the 
fact that the dynamics of lattice models of such systems 
generates fluctuations that can wash-out spatio-t emp oral 
patterns predicted by coarse-grained approaches [37H39j . 

In our model two species, which differ in the metabolic 
and reproduction rates, occupy lattice sites and consume 
the same renewable resources. In the context of tumor- 
growth problem, faster evolving species can be considered 
as tumor or cancer (c), the slower one as normal cells (n), 
and resources (p) become nutrients provided with blood. 
Each site of a d-dimensional lattice of linear size L can be 
empty or occupied by a nutrient, by a normal cell, or by a 
tumor cell (and we have implemented periodic boundary 
conditions). Moreover, we implement the intraspecific 
exclusion rule: no more than one cell of a given type can 
occupy a given site. It means that in our model each site 
can be in one of the eight states (p, n, c), where p, n, c = 
or 1. For example, (0, 1, 1) stands for the absence of a 
nutrient cell and the presence of both normal and tumor 
cells. The competition between normal and tumor cells 
for nutrients is modeled in a way that resembles prey- 
predator systems with n and c being two predator species, 
and p playing the role of preys. As a matter of fact, 
our model can be regarded as an extension of the prey- 



predator model that was already examined in the context 
of oscillatory and critical behaviour that such systems are 
known to exhibit [28|, H(| . 

The detailed rules of the model are specified below: 

• Choose a site randomly. 

• With probability r p update the site provided that 
there is a nutrient cell on this site (p = 1). In 
this case an update means an attempt to breed: 
one of the nearest neighbours of the selected site, 
not containing a nutrient cell, is chosen and a new 
nutrient cell is placed there. 

• With probability r„ update the site provided that 
there is a normal cell there (n = 1). In this case the 
normal cell survives only if there is also a nutrient 
cell on this site. If so, the normal cell consumes the 
nutrient cell {(p = 1) — > (p — 0)) and it attempts to 
breed (in a way analogous to that described above 
for a nutrient cell). 

• With probability r c update the site provided that 
there is a tumor cell there (c = 1). The tumor cell 
survives only if there is also a nutrient cell on this 
site. If so, the tumor cell consumes the nutrient 
cell and attempts to breed in an analogous way as 
a nutrient cell does. 

The update probabilities are the only control parame- 
ters of the model and they satisfy the obvious condition: 

r p + r n + r c = 1. (1) 

We assume that tumor cells have faster metabolism 
and reproduction rates than normal cells (i.e., r c > r n ). 
With Eq. ([1} and a fixed ratio of r c /r n , the behaviour of 
the model depends only on r p . 

To examine the model, we introduce eight densities 
Xpnc, which measure the average concentrations of the 
respective states in a system: 

x P nc = jd^2S(i,(p,n,c)), (2) 

i 

where summation is over all lattice sites i, and 
S(i, (p, n,c)) = 1 when the site i is in the state (p,n,c), 
and otherwise. Using the densities, we can calculate the 
concentrations of nutrient (x p ), normal (x n ), and tumor 
cells (x c ) as follows: 

Xp = XlOQ + Xn + X W i + Xm, 
x n = 2*010 + 2*110 + ^011 + 2*111, 
2*c = XqoI + ^lOl + ^011 + 2*1 11- (3) 

The densities x p , x n , and x c are quantities of our main 
interest. 

Although different metabolism and reproduction rates 
inclined us to interpret two species in our model as nor- 
mal and tumor cells, such an analogy must be taken with 
considerable care. Indeed, a number of features of real 
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normal-tumor competition is not properly reflected in 
our model. For example, normal and tumor cells can- 
not occupy the same place, proliferation of normal cells 
is extremely small comparing to the cancer cells, and the 
latter ones in some cases might even diffuse. Moreover, 
nutrients are much smaller than normal or tumor cells 
and they do not divide. We are aware of such shortcom- 
ings of our model, but the objective of understanding 
of competing systems from the perspective of statisitcal 
mechanics implied that a number of factors had to be 
omitted. As we will show in Section V some of the prop- 
erties of our model seem to be relatively robust, and there 
is a hope that more realistic systems will exhibit similar 
behaviour. 



III. MEAN-FIELD APPROXIMATION 

For the model defined by the above dynamical rules, 
one can write a Master equation and solve it using some 
decoupling procedure [H, [13] . The first step in this ap- 
proach is to sum over states of all but one site. As a 
result we obtain equations describing the time evolution 
of probabilities of a single-site configuration. Such equa- 
tions will contain probabilities of more complex configu- 
rations (i.e. two-site configurations) that need to be sub- 
sequently factorized (i.e., approximated with some prod- 
ucts of probabilities of single-site configurations). Even- 
tually, each term contributing to the evolution of single 
site probabilities will contain probability that a chosen 
site is in an appropriate state (for a given process) and 
the rate of that process (r p , r n or r c ). In some processes 
neighbouring sites also give some contributions. For our 
model one arrives at the following equations, describing 
the time evolution of the densities x. 
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(4) 



where the unit of time corresponds to a single (on aver- 
age) update per lattice site and 



x n = xi 10 + xm, 
x c = x W i + xm. 



(5) 



Moreover, the function f(x) denotes the average number 
of sites that had chosen to breed (with respect to x) at a 
given site 
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where w denotes the number of nearest neighbors (on 
cartesian lattices w = 2d). After simple algebra one can 
show that f(x) — . The quantities x n (x c ) denote 
the densities of sites with both normal (tumor) and nutri- 

obeys a similar equation, 



dxn 



dt 



ent cells. The derivative 
but it is simpler to use the obvious normalization condi- 
tion 
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i. 



(7) 



p,n,c— 0,1 



That mean-field equations @ are indeed consequences 
of dynamical rules can be seen also from less formal ar- 
guments. For example, the first term in the first equa- 
tion describes the increase of Xiqq due to the breeding 
of nutrient cells whose neighbour chosen for the place- 
ment of a newborn nutrient cell is empty. The factor 
f{x p )xooo = (I - x™ ) • T^aa. gives the probability that 
at least one of the neighbouring sites does not contain a 
nutrient cell ((I — x p )) and that the chosen neighbouring 
site (that does not contain a nutrient cell) is in the state 
(0, 0, 0) (the conditional probability of this equals ). 

We solved numerically the equations (HI and calculated 
the densities ([3]) . The results of our calculations together 
with those of Monte Carlo simulations are discussed be- 
low. Let us also notice that for some biologically-inspired 
lattice models (i.e., wound healing problem) other ap- 
proximate descriptions are possible based for example on 
Cahn-Hilliard equation [35]. 



IV. RESULTS OF MONTE CARLO 
SIMULATIONS 

Monte Carlo simulations are often used to study lattice 
models. Starting for each r p from a new random configu- 
ration (unless specified otherwise, we used configurations 
where each site with probability 1/4 was either empty or 
occupied by a nutrient, normal, or tumor cell), we mea- 
sured the steady-state densities ^ as functions of r p , for 
the fixed ratio r c /r n . Of course, before measurements we 
relaxed the system and monitored some densities, until it 
reached the steady state. To examine critical behaviour 
and the nature of some phase transitions in the model, 
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FIG. 1. Steady-state densities of nutrient (x p ), normal (x n ), 
and cancer (x c ) cells as functions of r p for d = 1 and 
r c /r n — 3. The simulations were made for L — 10 5 and 
simulation time t = 10 4 . Close to critical points, longer sim- 
ulations were performed with t = 10 6 . For each value of r p , 
an initial configuration was random and contained individu- 
als of all species. Solutions of the mean-field approximation 
(Q are also shown, but, as one can see, the agreement with 
simulations is rather poor. 



FIG. 2. The time dependence of the density of normal cells x n 
calculated for d = 1, r c /r n = 3, and (from top) r p — 0.205, 
0.2, 0.197, 0.195, 0.194, 0.193, and 0.19. The dotted line 
has a slope corresponding to the (1+1) directed percolation 
value, Sdp = 0.1595 41]. The inset shows the behaviour of 
x n in the vicinity of critical points located at r p — 0.1942 
(+) and r p — 0.6367 (x). In both cases, power-law behaviour 
is in a good agreement with the directed percolation value, 
I3 DP = 0.2765. 



we also measured the time dependence of some of these 
densities. Additional details concerning our simulations 
are provided in the following subsections. 



A. d=l 

The plot of the steady-state densities x p , x n , and x c as 
functions of r p for r c /r n = 3 is shown in Fig.[TJ For large 
r p mainly nutrient cells are updated, i.e., they breed if 
there is a place nearby. As a result, nutrients are easily 
available and both normal and tumor cells can coexist. 
For smaller values of r p , the shortage of nutrients affects 
mainly normal cells and for r p ~ 0.6367 they become 
extinct. For r p even smaller, only tumor cells can exist 
- although an initial configuration always contains indi- 
viduals of both species, normal cells quickly die out, and 
only tumor cells survive. Then, surprisingly, a drastic 
change takes place: for r p < 0.5195 these are normal 
cells that survive and tumor cells that become extinct. 
When r p is too small (for r p < 0.1942), nutrients do not 
reproduce sufficiently fast and both normal and tumor 
cells become extinct (and only nutrients survive). 

As might be expected, for low-dimensional systems the 
solution predicted by the mean-field approximation Q 
with w = 2 considerably differs from Monte Carlo simu- 
lations (see Fig. [1]). 

The limits of existence of normal or tumor cells are 
marked by phase transitions in the model. To exam- 
ine their nature, we calculated the time-dependence of 



the average density of species that becomes extinct at 
the transition point (usually the averages are taken over 
10 2 independent runs, which start from different random 
initial configurations). One expects that at the critical 
point these quantities show power-law decay x(t) ~ t~ s , 
where the characteristic exponent 5 exhibits some degree 
of universality [4l| . The results of our calculations for the 
transition at r p — 0.1942 are shown in Fig. [21 and from 
the behaviour at the critical point we estimate 5 ~ 0.16. 
This value is very close to the corresponding exponent 
in the directed percolation (DP) problem, 8dp = 0.1595 
[4lj . An additional argument that the model belongs to 
the DP universality class comes from the analysis of the 
steady-state density of normal cells x n close to the critical 
point. We find that x n ~ (r p — 0.1942) 19 , where /3 is very 
close to the directed percolation value, Pop — 0.2765 [41[ 
(see the inset in Fig. The transition at r p — 0.1942 
is similar to some other transitions in models with ab- 
sorbing states [29L l4lj . The critical behaviour of the DP 
universality class is typical for a large class of models with 
a single absorbing state and our model also falls into this 
class [42l ]. 

Similar calculations were performed for the transition 
at r p = 0.6367, where normal cells become extinct. We 
show only the steady-state values of the density of normal 
cells (inset in Fig. [5]) and they suggest that also in this 
case the model most likely belongs to the DP universal- 
ity class. The time dependence of x n at the critical point 
(not presented) yields the estimation of 6, which is also 
consistent with the DP value. Let us notice, however, 
an important difference between the phase transition at 
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FIG. 3. The time dependence of the density of tumor cells 
x c calculated for d = 1, r c /r n — 3, and (from top) r p = 0.54, 
0.522, 0.52, 0.5195, 0.519, 0.518, 0.517, and 0.515. 

0.1942, where the model enters an absorbing state, and 
the one at 0.6367, where instead of entering the absorbing 
state, the system remains in the active phase with tumor 
and nutrient cells (but without normal cells) . This phase 
transition provides yet another example that the DP crit- 
ical behaviour appears for a wider class of models than 
those with a single absorbing state. However, a complete 
understanding of this issue is not yet achieved |43l - l45| . 

A different behaviour is observed for the transition 
at r p — 0.5195. The steady-state densities exhibit pro- 
nounced jumps (Fig.UJ and that suggests a discontinuous 
transition at this point. The plot of time dependence of 
the tumor cells confirms such behaviour (Fig. [3]). Indeed, 
no power-law decay of x c as a function of time is seen at 
the critical point. 

The above-described behaviour seems to be generic in 
one-dimensional version of our model. We observed the 
same sequence of phase transitions (DP-discontinuous- 
DP), albeit at other values of r p , for r c /r n ranging from 
1.1 to 10. It would be interesting to find a reason of such 
a robust behaviour and we will return to this problem in 
the next subsection. 

In studies on tumor growth, it is often important to 
examine the time evolution of frequency or distribution 
of tumor cells Q. The observed patterns of invasion of 
tumor cells sometimes resemble a scenario known to ecol- 
ogists as "seed and soil", where a new species colonizes 
a new habitat [46]. In our model a qualitatively similar 
analysis can be made. In Fig. U we show an example 
of tumor remission. In the initial configuration, 80% of 
the central sites were occupied by nutrient and tumor 
cells (1,0,1) and the rest by nutrient and normal cells 
(1,1,0). The simulations were performed for r c /r n = 3 
and r p = 0.48, for which we have only nutrient and nor- 
mal cells in the steady state. In Fig. @]a clear extinction 
of tumor cells takes place. An opposite effect is seen in 
Fig. [5] where the initial configuration contained only a 
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FIG. 4. Tumor remission. The time evolution of distribution 
of tumor (top) and normal (bottom) cells calculated for d = 1, 
r c /r n = 3, and r v = 0.48. Sites occupied by tumor (top) or 
normal (bottom) cells are plotted with black dots. 

small cluster of 10 tumor cells, while r v = 0.55. Under 
such conditions a tumor recurrence can be observed with 
normal cells becoming gradually extinct. 

B. d=2 

We made similar calculations for the two-dimensional 
model. For r c /r n — 3, steady-state densities are shown 
in Fig. |51 Although there are similarities to the d = 
1 case, there are also some qualitative differences. In 
particular, a single discontinuous transition that reverses 
the extinction of tumor and normal cells is replaced in 
the d = 2 version by two transitions and a narrow range 
of a coexistence: 0.2947 < r p < 0.3441 (see Fig. [6]). A 
more detailed analysis of the densities close to the phase 
transitions shows that three of these phase transitions are 
likely to belong to the DP(2+1) universality class (see 
inset in Fig. However, the scaling of the density of 
normal cells close to the phase transition at r p — 0.3441 
seems to scale with the exponent (3 ~ 0.32, which is much 
smaller than the DP(2+1) exponent f3 DP = 0.584 f4ljj . 

An additional argument that this phase transition 
might be described by exponents different than DP(2+1) 
is presented in Fig. [7j which shows on the log-log scale 
the time dependence of the density of normal cells x n . At 
the transition point (thick line) , the asymptotic decay of 
x n seems to be characterized by the exponent S w 0.25, 
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FIG. 7. The time dependence of the density of normal cells 
x n calculated for d — 2, r c /r n = 3, and (from top) r p = 
0.341, 0.343, 0.3435, 0.3437, 0.3439, 0.344. 0.3441 (thick line), 
0.3442, 0.3445, 0.345, 0.347, and 0.35. The dotted line has a 
slope corresponding to the (2+1) directed percolation value, 
S D p = 0.451 0. 



FIG. 5. Tumor recurrence. The time evolution of the distri- 
bution of tumor (top) and normal (bottom) cells calculated 
for d = 1, r c /r n = 3, and r v — 0.55. Sites occupied by tumor 
(top) or normal (bottom) cells are plotted with black dots. 
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FIG. 6. Steady-state densities of nutrient (%,), normal (x n ), 
and cancer (x c ) cells as functions of r p for d — 2 and r c /r n — 
3. Simulations were made for L — 500 and simulation time 
t — 10 4 . Close to critical points, longer simulations were 
performed with t — 10 6 . The inset shows the behaviour of 
x n and x c in the vicinity of critical points located at (from 
top) r p = 0.3441, 0.4249, 0.2947, and 0.0299. Except for r p = 
0.3441, the power-law behaviour is in a good agreement with 
the (2+1) directed percolation value, Pdp = 0.584 (dashed 
line) HJ]. 



which is different from 8 dp = 0.451 [4l|. We do not 
present numerical data here, however, for the remaining 
three transitions the estimated exponent 5 was very close 
to the DP(2+1) value. 

Despite such a discrepancy with the expected DP(2+1) 
universality class, our estimations of critical exponents 
must be taken with some care. Namely, it is known that 
such estimations for nonequilibrium phase transitions are 
often very difficult and we cannot exclude that much 
more extensive numerical calculations will modify our es- 
timations. Let us also notice that the phase transition at 
r p = 0.3441 has an active subcritical phase (nutrient and 
tumor cells) , and such a feature might perhaps cause nu- 
merical difficulties or even change the critical behaviour 
(however, the phase transition at r p = 0.4249 seems to 
have the DP(2+1) critical exponents, but also have an 
active subcritical phase). 

Tumors in two-dimensional models have much richer 
morphology than those in the one-dimensional case. 
Figs. [5] El show the distribution of tumor cells that grew 
from a small tumor seed surrounded by normal and nu- 
trient cells. Let us notice that for r p = 0.34, which is rel- 
atively far from the limit of existence of tumor (0.2947), 
the tumor remains compact and dense (Fig. [H]). For 
r p = 0.304, i.e., very close to such a limit of existence, the 
shape of tumor is much more irregular (Fig. [8]) . Morphol- 
ogy of tumors is very important from the clinical point of 
view, however, a precise comparison with existing data is 
very often difficult 47j. Nevertheless, certain morpholo- 
gies can be associated with some forms of cancer [16| . 

Let us notice that the coexistence of tumor and nor- 
mal cells in the range 0.2947 < r p < 0.3441 is very fragile 
since this interval is rather narrow and both species re- 
main close to their limits of existence. To have some in- 
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FIG. 8. The snapshot distribution of tumor cells calculated 
for d — 2, r c /r n = 3, r p = 0.304 and simulation time t — 10 4 . 
The initial configuration contained a small seed of tumor cells 
surrounded by normal and nutrient cells. 

sight into the nature of this coexistence, we present the 
distribution of nutrient and tumor cells (Fig. [TU|). We 
can observe that tumor cells seem to concentrate close to 
empty regions. Apparently, under such conditions, in the 
interior they lose the competition with normal cells. It 
is tempting to speculate that such structures might have 
a rather short lifetime in one-dimensional systems and 
that is why there is no coexistence for the d = 1 model 
except in the large r p regime (see Fig. [I]). 

It turns out, however, that the behaviour in the d = 2 
version, which we have shown in Fig. [51 is not entirely 
generic and to some extent depends on the ratio r c /r n . 
Namely, numerical simulations (which we do not present) 
show that for r c /r n = 1.5 and 1.1 some of the continuous 
transitions become discontinuous, as in one-dimensional 
system, but with a narrow range of r p for which both tu- 
mor and normal cells coexist. For larger r c /r n (we made 
simulations for r c /r n — 5) all phase transitions remain 
continuous. Taking into account appearance of discon- 
tinuous transitions one cannot exclude that the phase 
transition for r c /r n = 3 at r p = 0.3441 is actually weakly 
discontinuous and that is why our estimations of critical 
exponents do not agree with the DP universality class 
exponents. 

C. d=3 

We made numerical simulations also for the d = 3 
version. For r c /r n — 3, the steady-state densities as 
functions of r p are shown in Fig. 111! We did not esti- 




FIG. 9. The snapshot distribution of tumor cells calculated 
for d = 2, r c /r n = 3, r v = 0.34, and simulation time t = 2-10 3 . 
The initial configuration contained a small seed of tumor cells 
surrounded by normal and nutrient cells. 




FIG. 10. The snapshot distribution of nutrient (•) and tumor 
(*) cells calculated for d — 2, r c /r n — 3, and r v — 0.305. 
Tumor cells seem to concentrate close to empty regions - ap- 
parently in the interior they lose the competition with normal 
cells. 
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FIG. 11. Steady-state densities of nutrient (%,), normal (x„), 
and cancer (x c ) cells as functions of r p for d — 3 and r c /r n = 
1.5. Simulations were made for L = 50 and simulation time 
i = 10 4 . Close to critical points, longer simulations were 
performed with L = 80 and t = 10 5 . The solutions of the 
mean-field approximation ([4]) are also shown. For 0.05 < r p > 
0.22 normal cells and nutrients show oscillatory behaviour and 
in such a case the plotted values should be interpreted as time 
averages rather than steady-state values. 



mate the critical exponents, but by analogy with low- 
dimensional versions, we expect that continuous transi- 
tions belong to the DP (3+1) universality class. Let us 
also notice a good agreement with numerical solutions of 
the mean-field approximation Q with w = 6. With a 
pronounced discontinuous transition seen in Fig. 111! the 
three-dimensional model resembles the one-dimensional 
version. It seems to us, however, that there are some im- 
portant differences between these two cases. Namely, we 
noticed that for d = 3 the location of the discontinuous 
transition depends on the concentration of species in the 
initial configuration. For example, for initial configura- 
tions containing much more tumor than normal cells, the 
location of this transition is shifted toward much smaller 
r p (see Fig.fl2|). The dependence on the initial configu- 
ration appears only in some vicinity of the discontinuous 
transition. In Fig. Q2J for r p < 0.18 or 0.33 < r p < 1 , 
the steady-state densities are within numerical error the 
same as those in Fig. [TT] Such a dependence on the ini- 
tial configuration is also reproduced in the solutions of 
the mean-field equations (j4j. In our opinion it is likely 
that these are not only densities of species but also their 
spatial distribution that determines location of disconti- 
nuities. More detailed analysis of such effects would be 
certainly desirable. 

On the other hand, numerical simulations (not pre- 
sented) show that in one-dimensional systems location of 
discontinuity is essentially independent on the initial con- 
figuration (provided that it contains a non-zero concen- 
tration of normal and tumor cells, and nutrients). The 
simulations show also that in the two-dimensional ver- 



FIG. 12. Steady-state densities of nutrient (x p ), normal (x„), 
and cancer (x c ) cells as functions of r v for d = 3 and r c /r n = 
1.5. Simulations were made for L = 50 and simulation time 
t — 10 4 . In the initial configuration there was much more 
tumor than normal cells (x p = 0.2, x„ — 0.05, and x c — 0.65). 
For r v < 0.18 or 0.33 < r p < 1, the steady-state densities are 
within numerical error the same as those in Fig. [TTJ 

sion, for r c /r n = 1.5 and 1.1, location of discontinuity 
depends on the initial configuration. For r c /r n = 3 or 5 
all transitions seem to be continuous and no dependence 
on the initial configuration was detected. It would be 
desirable to know whether there are some general physi- 
cal arguments about the role of dimension in determining 
such a behaviour of our model. 



D. d>3 

Systems with long-range interactions or fast mixing 
due to diffusion are often well approximated by infinite- 
dimensional systems. It is thus of interest to examine 
models with d > 3. Since Monte Carlo simulations are 
very demanding for such models, we performed only cal- 
culations within the mean-field approximation (|TJ). This 
approximation was quite satisfactory already for d = 3 
(Fig. [TTJ , and we expect that for d > 3 the solutions of 
will be even more accurate. Our calculations show 
that for d = 4, 5, . . . , 10 the steady-state densities behave 
similarly to the d = 3 version. Location of the disconti- 
nuity was also found to depend on the initial condition. 

One of the interesting problems in lattice models con- 
cerns the existence of oscillations and numerous examples 
ranging from hetero- or homogeneous catalysis to various 
ecological systems have already been examined with this 
respect Q • On general grounds one expects that fluctu- 
ations in low-dimensional systems are strong enough to 
destroy such temporal oscillations in the thermodynamic 
limit. An important question is: what is the critical di- 
mension d c above which such oscillations will survive in 
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FIG. 13. Boundary of oscillatory regime calculated using 
the mean-field approximation @ for r c /r n = 3. Oscillations 
occur only in the tumor-free phase of the model (x c — 0). 

the thermodynamic limit? Grinstein et al. [49| presented 
certain arguments indicating that d c = 2. Numerical 
results for both synchronous jsjj [EH and asynchronous 
models [28| seem to confirm that temporal oscillations 
exist only for d>2. 

The mean-field approximation (j4j in some cases also 
predicts that solutions are oscillatory. For a given d, we 
analysed the range of r p where the oscillatory behaviour 
was seen, and the corresponding plot is shown in Fig. 1131 
Using this approximation, we estimate that d c ~ 1.8 and 
that result is very close to the already mentioned predic- 
tion d c = 2. 

V. RANDOMLY SUPPLIED NUTRIENTS 

As we already mentioned, nutrient cells in our model 
evolve according to rules typical rather to some living 
species. In particular, they divide and place their off- 
springs on some empty lattice sites. Such rules might 
be justified in modeling multi-species ecosystems where 
nutrients could be considered as prey-type species. How- 
ever, in the context of tumor growth problems division 
of nutrient cells is not realistic and different rules should 
be used. In the present section we examine a modifica- 
tion of our model, where nutrient cells are supplied ran- 
domly. More precisely, we left unchanged all the rules 
of the model (sect. II) except the second rule that now 
becomes 

• With probability r p place a nutrient cell at the cho- 
sen site, provided that there is no such a cell there. 

Simulations of such a model for d = 1,2 and 3 show 
that this modification has only a minor effect on the be- 
haviour of the model (Fig 03}. Indeed, the regions of 
existence of and coexistence resemble those in the previ- 
ous d = 3 version (Fig. ITT]) . We do not present our data 



FIG. 14. Steady-state densities of nutrient (x p ), normal 
(x n ), and cancer (x c ) cells as functions of r v for d = 3 and 
r a /r n — 1.5. Simulations were made for the version with ran- 
dom supply of nutrient cells and with L = 50 and simulation 
time t = 10 4 . 



but simulations for d = 1 and 2 also give qualitatively 
very similar results. 



VI. CONCLUSIONS 

In the present paper we analysed a simple lattice 
model of two species competing for common renewable 
resources. Although various models have already been 
examined in the literature, little attention was paid to 
general properties of such systems as seen from the many- 
body perspective. The objective of the present paper was 
to examine a competition of species from this very per- 
spective of statistical mechanics. The model might de- 
scribe the competition of two predator species for a prey 
but to some extent also of tumor and normal cells for 
some nutrients. It is the latter interpretation that was 
mainly used in our paper. 

The introduced model shows rich behaviour. It ex- 
hibits continuous transitions, which most likely belong 
to the directed percolation universality class. In some 
cases, however, the identification of the universality class 
was inconclusive and it was suggested that a novel criti- 
cal behaviour might appear in this model. There are also 
discontinuous transitions in the system. In the d = 1 ver- 
sion the extinction of tumor cells is always discontinuous, 
but for d > 1 versions both continuous and discontinuous 
extinctions of tumor cells appear. It would be desirable to 
explain such dimension-dependent behaviour using some 
more general physical arguments. 

Contrary to the Cause exclusion principle, our model 
show that in some cases coexistence of predator species 
that depend on the same prey is possible. In one- 
dimensional systems coexistence appears only when prey 
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is sufficiently abundant but higher dimension facilitates 
coexistence. It would be interesting however to general- 
ize such observations examining for example models with 
interspecific exclusion. 

Especially in the context of tumor growth problem it 
would be desirable to examine in more details the role 
of the initial state. Preliminary results show that for 
d > 1 the final state depends on the initial concentrations 
of normal and tumor cells, but most likely their spatial 
distribution plays a role as well. Finding distributions or 
configurations that can suppress spreading of tumor cells 
would be particularly valuable. 

In the present model, species do not change their be- 
haviour in time. A more comprehensive description of 



ecosystems as well as tumors [3J, in addition to ecolog- 
ical, should take into account also evolutionary aspects. 
Analysis of the evolutionary, multi-species extension of 
the present model might further contribute to a better 
understanding of such complex, but immensely impor- 
tant systems. 
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